#ifndef EQUATIONSOLVER
#define EQUATIONSOLVER

#include "Function.hpp"
#include <cstdio>
#include <math.h>

// define sgn()
inline int sgn(double x){
    return (x>0) - (x<0);
}

class EquationSolver{
protected:
    const Function & F;
public:
    EquationSolver(const Function& F) : F(F) {}
    virtual double solve() = 0;
};

class Bisection_Method : public EquationSolver {
private:
    double a, b;
    double eps, delta;
    int Maxiter;
public:
    Bisection_Method(const Function &F, double a, double b, 
        double eps = 1e-7, double delta = 1e-6, int Maxiter = 50) :
        EquationSolver(F), a(a), b(b), eps(eps), delta(delta), Maxiter(Maxiter) {}
    
    virtual double solve() {
        double h = b - a;
        if (h <= 0) {
            printf("Input error: b should be bigger than a.\n");
            return NAN; // indicate error
        }

        // u,v are undefined
        double u = F(a), v = F(b);
        if(isinf(u)) {
            u = F(a + 0.0001);
        }
        if(isinf(v)) {
            v = F(b - 0.0001);
        }
        
        if(sgn(u)==sgn(v)) {
            printf("Input error: sgn(f(a)) should be different from sgn(f(b)).\n");
            return NAN; // indicate error
        }

        double c;
        for (int k = 0; k < Maxiter; k++)
        {
            h = h/2;
            c = a + h;

            if (fabs(h) < delta) break;

            double w = F(c);
            if (fabs(w) < eps) break;

            if (sgn(w) == sgn(u)){
                a = c; 
            } 
        }
        return c;   // Return the approximate root.
    }
};

class Newton_Method : public EquationSolver {
private:
    double x0;
    double eps;
    int Maxiter;
public:
    Newton_Method(const Function &F, double x0, 
        double eps = 1e-7, int Maxiter = 8) :
        EquationSolver(F), x0(x0),  eps(eps), Maxiter(Maxiter){}
    
    virtual double solve() {
        double x = x0;
        double u;
        for (int k = 0; k < Maxiter; k++)
        {
            u = F(x);
            if(fabs(u) < eps) break;
            x = x - u/F.derivative(x);
        }
        return x;   // Return the approximate root.
    }
};

class Secant_Method : public EquationSolver {
private:
    double x0, x1;
    double eps, delta;
    int Maxiter;
public:
    Secant_Method(const Function &F, double x0, double x1,
        double eps = 1e-7, double delta = 1e-6, int Maxiter = 20) :
        EquationSolver(F), x0(x0), x1(x1), eps(eps), delta(delta), Maxiter(Maxiter) {}
    
    virtual double solve() {
        double x_now = x1, x_last = x0;
        double u = F(x_now), v = F(x_last);
        double s;
        for (int k = 2; k < Maxiter; k++)
        {
            s = (x_now - x_last)/(u-v);
            x_last = x_now;
            x_now = x_now - u*s;
            v = u;
            if(fabs(x_now-x_last) < delta) break;
            u = F(x_now);
            if(fabs(u) < eps) break;
        }
        return x_now;   // Return the approximate root.
    }
};

#endif